Chapter 12 

Introduction to Simulation Using 

MATLAB 

A. Rakhshan and H. Pishro-Nik 


12.1 Analysis versus Computer Simulation 

A computer simulation is a computer program which attempts to represent the real world based 
on a model. The accuracy of the simulation depends on the precision of the model. Suppose that 
the probability of heads in a coin toss experiment is unknown. We can perform the experiment 
of tossing the coin n times repetitively to approximate the probability of heads. 

Number of times heads observed 
Number of times the experiment executed 

However, for many practical problems it is not possible to determine the probabilities by exe¬ 
cuting experiments a large number of times. With today’s computers processing capabilities, 
we only need a high-level language, such as MATLAB, which can generate random numbers, to 
deal with these problems. 

In this chapter, we present basic methods of generating random variables and simulate prob¬ 
abilistic systems. The provided algorithms are general and can be implemented in any computer 
language. However, to have concrete examples, we provide the actual codes in MATLAB. If you 
are unfamiliar with MATLAB, you should still be able to understand the algorithms. 

12.2 Introduction: What is MATLAB? 

MATLAB is a high-level language that helps engineers and scientists find solutions for given 
problems with fewer lines of codes than traditional programming languages, such as C/C++ 
or Java, by utilizing built-in math functions. You can use MATLAB for many applications 
including signal processing and communications, finance, and biology. Arrays are the basic data 
structure in MATLAB. Therefore, a basic knowledge of linear algebra is useful to use MATLAB 
in an effective way. Here we assume you are familiar with basic commands of MATLAB. We 
can use the built-in commands to generate probability distributions in MATLAB, but in this 
chapter we will also learn how to generate these distributions from the uniform distribution. 
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12.3 Discrete and Continuous Random Number Generators 

Most of the programming languages can deliver samples from the uniform distribution to us 
(In reality, the given values are pseudo-random instead of being completely random.) The rest 
of this section shows how to convert uniform random variables to any other desired random 
variable. The MATLAB code for generating uniform random variables is: 


U = rand; 


which returns a pseudorandom value drawn from the standard uniform distribution on the open 
interval (0,1). Also, 


U = random, n); 


returns an rn-by-n matrix containing independent pseudorandom values drawn from the standard 
uniform distribution on the open interval (0,1). 


12.3.1 Generating Discrete Probability Distributions from Uniform Distri¬ 
bution 


Let’s see a few examples of generating certain simple distributions: 


Example 1. (Bernoulli) Simulate tossing a coin with probability of heads p. 

Solution: Let U be a Uniform(0,l) random variable. We can write Bernoulli random variable 
X as: 


j 1 U < p 
\ 0 U > p 


Thus, 


P{H) = P(X = 1) 
= P(U < p) 
= P 


Therefore, X has Bernoulli(p ) distribution. The MATLAB code for Bernoulli^ 0.5) is: 


p = 0.5; 

U = rand; 

X = (U < p); 


Since the “rand” command returns a number between 0 and 1, we divided the interval [0,1] 
into two parts, p and 1 — p in length. Then, the value of X is determined based on where the 
number generated from uniform distribution fell. 
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Example 2. (Coin Toss Simulation) Write codes to simulate tossing a fair coin to see how the 
law of large numbers works. 

Solution: You can write: 


n = 1000; 

U = rand(l, n); 

toss = (U < 0.5); 

a = zeros(n + 1); 

avg = zeros{n ); 

for i = 2 : n + 1 

a(i) = a{i — 1) + toss{i — 1); 

avg(i - 1) = a(i)/(i - 1); 

end 

plot(avg) 


If you run the above codes to compute the proportion of ones in the variable “toss,” the result 
will look like Figure 12.5. You can also assume the coin is unbiased with probability of heads 
equal to 0.6 by replacing the third line of the previous code with: 


toss = (U < 0.6); 
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Figure 12.1: MATLAB coin toss simualtion 


Example 3. (Binomial) Generate a Binomial (50,0.2) random variable. 
Solution: To solve this problem, we can use the following lemma: 





4 CHAPTER 12. INTRODUCTION TO SIMULATION USING MATLABA. RAKHSHAN AND H. PISHRO-NIK 


Lemma 1. If X\,X 2 , ..., X n are independent Bernoulli(p) random variables, then the random 
variable X defined by X = X\ + X 2 + ... + X n has a Binomial{n,p) distribution. 

To generate a random variable X ~ Binomial(n,p), we can toss a coin n times and count the 
number of heads. Counting the number of heads is exactly the same as finding X\+X 2 + ...+X n , 
where each X{ is equal to one if the corresponding coin toss results in heads and zero otherwise. 

Since we know how to generate Bernoulli random variables, we can generate a Binomial{n,p) 
by adding n independent Bernoulli(p ) random variables. 


P = 0-2; 
n = 50; 

U = rand(n , 1); 

X = sum(U < p); 


Generating Arbitrary Discrete Distributions 


In general, we can generate any discrete random variables similar to the above examples using 
the following algorithm. Suppose we would like to simulate the discrete random variable X with 
range R\ = {xi,x 2 , ...,x n } and P(X = Xj) = pj , so YljPj = 1- 

To achieve this, first we generate a random number U (i.e., U ~ Uniform^ 0,1)). Next, we 
divide the interval [0,1] into subintervals such that the jth subinterval has length pj (Figure 
12.2). Assume 


X=l 


Xo 

if 

0 

a, 

V 

b 

X\ 

if 

(po <U < po+Pi) 

Xj 

if 

(ESb < V < Ei 


In other words 

X = Xj if F{xj- 1 ) <U< F(xj), 
where F{x) is the desired CDF. We have 


P(X = Xj ) = P 


3 -1 3 \ 

^2pk <u < 

.k =0 k =0 / 


= Pi 


|P0^1 | P2 | P3 | ••• Pj \ 
0 1 


Figure 12.2: Generating discrete random variables 
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Example 4. Give an algorithm to simulate the value of a random variable X such that 

P(X = 1) = 0.35 
P(X = 2) = 0.15 
P(X = 3) = 0.4 
P(X = 4) = 0.1 

Solution: We divide the interval [0,1] into subintervals as follows: 

A 0 = [0,0.35) 

Ai = [0.35,0.5) 

A 2 = [0.5, 0.9) 

A 3 = [0.9,1) 

Subinterval Aj has length p^. We obtain a uniform number U. If U belongs to A.;, then X = x*. 

P(X = Xi ) = P(U g M) 

= Pi 


P = [0.35,0.5,0.9,1]; 

X = [1,2,3,4]; 

counter = 1; 
r = rand ; 

while(r > P(counter)) 
counter = counter + 1; 
end 

X ( counter ) 


12.3.2 Generating Continuous Probability Distributions from the Uniform 
Distribution- Inverse Transformation Method 

At least in principle, there is a way to convert a uniform distribution to any other distribution. 
Let’s see how we can do this. Let U ~ Uniform^ 0,1) and F be a CDF. Also, assume F is 
continuous and strictly increasing as a function. 

Theorem 1. Let U ~ Uniform^ 0,1) and F be a CDF which is strictly increasing. Also, 
consider a random variable X defined as 


X = F _1 (17). 


Then. 


A ~ F (The CDF of X is F) 




6 CHAPTER 12. INTRODUCTION TO SIMULATION USING MATLABA. RAKHSHANAND H. PISHRO-NIK 


Proof: 


P(X < x) = P(F~ 1 (U) < x ) 

= P(U < F(x)) (increasing function) 

= F(x) 

Now, let’s see some examples. Note that to generate any continuous random variable X with 
the continuous cdf F, F _1 ([7) has to be computed. 

Example 5. (Exponential) Generate an Exponential(l) random variable. 

Solution: To generate an Exponential random variable with parameter A = 1, we proceed 
as follows 


F(x) = 1 — e x x > 0 
U Uniform,(0, 1) 

X = F~ l (U) 

= - ln(l - U) 

X ~ F 


This formula can be simplified since 


1 — U ~ Uniform{ 0,1) 


U, l-U 


Figure 12.3: Symmetry of Uniform 


Hence we can simulate X using 


X = - In (U) 


U = rand ; 

X = -log(Uy, 


Example 6. (Gamma) Generate a Gamma(20,l) random variable. 

Solution: For this example, F -1 is even more complicated than the complicated gamma cdf 
F itself. Instead of inverting the CDF, we generate a gamma random variable as a sum of n 
independent exponential variables. 

Theorem 2. Let X\, X 2 , ■ ■ ■ , X n be independent random variables with X{ ~ Exponential (A). 
Define 


Y = X 1 + X 2 + ■ • • + X n 
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By the moment generating function method, you can show that Y has a gamma distribution 
with parameters n and A, i.e., Y ~ Gamm,a(n, A). 

Having this theorem in mind, we can write: 


n = 20; 
lambda = 1; 

X = (—1/lambda) * sum/log(rand(n , 1))); 


Example 7. (Poisson) Generate a Poisson random variable. Hint: In this example, use the fact 
that the number of events in the interval [0, t] has Poisson distribution when the elapsed times 
between the events are Exponential. 

Solution: We want to employ the definition of Poisson processes. Assume N represents the 
number of events (arrivals) in [0,t]. If the interarrival times are distributed exponentially (with 
parameter A) and independently, then the number of arrivals occurred in [0,t], N, has Poisson 
distribution with parameter A t (Figure 12.4). Therefore, to solve this problem, we can repeat 
generating Exponential(X) random variables while their sum is not larger than 1 (choosing 
t = 1). More specifically, we generate Exponential(X) random variables 


Tt = In (Ui) 


by first generating uniform random variables U/'s. Then we define 


X = max { j '■ T\ + ■ ■ ■ + T.j < 1} 


The algorithm can be simplified: 




max 


-1 

X 


ln(C/i • 


Uj)< 1 
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Lambda = 2; 
i = 0; 

U = rand ; 

Y = —(1/ Lambda) * log(JJ)\ 
sum = Y ; 

while(sum <= 1) 

U = rand; 

Y = —(1/ Lambda) * log(U ); 
sum = sum + Y ; 

i = i + 1; 
end 
X = i- 


I Exp(X) ! Exp(X) j Exp(X) ! Exp(X) ^ | 

0 1 
X= Maximum number of exponential random variables 

Figure 12.4: Poisson Random Variable 


To finish this section, let’s see how to convert uniform numbers to normal random variables. 
Normal distribution is extremely important in science because it is very commonly occuring. 

Theorem 3. (Box-Muller transformation) We can generate a pair of independent normal vari¬ 
ables (Z i, Z‘i) by transforming a pair of independent Uniform^ 0,1) random variables {U\, U 2 ) 
[!]• 


J Z\ = \J —2 In Vi cos(27t[/2) 

I Z 2 = \/—2 In U\ sin(27r[/2) 

Example 8. (Box-Muller) Generate 5000 pairs of normal random variables and plot both 
histograms. 

Solution: We display the pairs in Matrix form. 


r = rand{ 5000, 2); 

n = sqrt {—2 * log(r( 1))) * [1,1]. * [cos(2 * pi * r(:, 2)), sin (2 * pi * r(:, 2))]; 
hist(n) 
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1500 



Figure 12.5: Histogram of a pair of normal random variables generated by Box-Muller transfor¬ 
mation 

12.4 MATLAB Commands for Special Distributions 

In this section, we will see some useful commands for commonly employed distributions. To be 
as precise as possible, we repeat the description of the commands from MATLAB help [2]. 

12.4.1 Discrete Distributions 

- Binomial Distribution: 

Y = binopdf(X,N,P) 

computes the binomial pdf at each of the values in X (vector) using the corresponding 
number of trials in N and probability of success for each trial in P. Y, N, and P can be 
vectors, matrices, or multidimensional arrays that all have the same size. A scalar input 
is expanded to a constant array with the same dimensions of the other inputs. 


Y = binocdf(X,N,P) 

computes a binomial cdf at each of the values in X using the corresponding number of 
trials in N and probability of success for each trial in P. X, N, and P can be vectors, ma¬ 
trices, or multidimensional arrays that are all the same size. A scalar input is expanded 
to a constant array with the same dimensions of the other inputs. The values in N must 
all be positive integers, the values in X must lie on the interval [0,N], and the values in P 
must lie on the interval [0, 1]. 

R = binornd(N,P) 

generates random numbers from the binomial distribution with parameters specified by 
the number of trials, N, and probability of success for each trial, P. N and P can be vectors, 
matrices, or multidimensional arrays that have the same size, which is also the size of R. 
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A scalar input for N or P is expanded to a constant array with the same dimensions as 
the other input. 

- Poisson Distribution 

Y = poisspdf(X,lambda) 

computes the Poisson pdf at each of the values in X using mean parameters in lambda. 

P = poisscdf(X,lambda) 

computes the Poisson cdf at each of the values in X using the corresponding mean param¬ 
eters in lambda. 

R = poissrnd(lambda) 

generates random numbers from the Poisson distribution with mean parameter lambda. 

- Geometric Distribution 

Y = geopdf(X,P) 

computes the geometric pdf at each of the values in X using the corresponding probabili¬ 
ties in P. 

Y = geocdf(X,P) 

computes the geometric cdf at each of the values in X using the corresponding probabilities 
in P. 

R = geornd(P) 

generates geometric random numbers with probability parameter P. P can be a vector, a 
matrix, or a multidimensional array. 

12.4.2 Continuous Distributions 

- Normal Distribution: 

Y = normpdf(X, mu, sigma) 

computes the pdf at each of the values in X using the normal distribution with mean mu 
and standard deviation sigma. 


P = normcdf(X,mu,sigma) 

computes the normal cdf at each of the values in X using the corresponding mean mu and 
standard deviation sigma. 
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nlogL = normlike(params,data) 

returns the negative of the normal log-likelihood function. 

R = normrnd(mu, sigma) 

R = randn 

generates random numbers from the normal distribution with mean parameter mu and 
standard deviation parameter sigma. 


- Exponential Distribution: 

Y = exppdf(X,mu) 

returns the pdf of the exponential distribution with mean parameter mu, evaluated at the 
values in X. 

P = expcdf(X,mu) 

computes the exponential cdf at each of the values in X using the corresponding mean 
parameter mu. 

R = exprnd(mu) 

generates random numbers from the exponential distribution with mean parameter mu. 


12.5 Exercises 

1. Write MATLAB programs to generate Geometric(p) and Negative Binomial(i,p) random 
variables. 

Solution: To generate a Geometric random variable, we run a loop of Bernoulli trials until 
the first success occurs. K counts the number of failures plus one success, which is equal 
to the total number of trials. 


K= 1 ; 

p = 0.2; 

while{rand > p) 
K = K + 1; 
end 
I< 


Now, we can generate Geometric random variable i times to obtain a Negative Binomial (i,p) 
variable as a sum of i independent Geometric (p) random variables. 
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2. (Poisson) Use the algorithm for generating discrete random variables to obtain a Poisson 
random variable with parameter A = 2. 

Solution: We know a Poisson random variable takes all nonnegative integer values with 
probabilities 

Pi = P(X = Xi) = e~ x - - for i = 0, 1,2, ••• 

i\ 

To generate a Poisson(X), first we generate a random number U. Next, we divide the 
interval [0,1] into subintervals such that the jth subinterval has length pj (Figure 12.2). 
Assume 

x 0 if (U < po ) 

xi if {po < U < po+Pi) 

X = 

x 3 if (£fc=JPfc < U < El= 0 Pk) 


Here x\ = i — 1, so 


X = i if po H- Pi-i < U < po-\ -h pi-i + Pi 

F{i - 1) < U < F{i) F is CDF 
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lambda = 2; 
i = 0; 

U = rand ; 

cdf = exp(—lambda); 
while(U >= cdf ) 
i = i + 1; 

cdf = cdf + exp(—lambda) * lambda A i/gamma(i + 1); 
end] 

X = i] 


3. Explain how to generate a random variable with the density 

f(x) = 2.5 xyfx for 0 < x < 1 

if your random number generator produces a Standard Uniform random variable U. Hint: 
use the inverse transformation method. 

Solution: 


F x (X) = xl=U (0<z<l) 

X = ui 


U = rand] 
X = 17§; 

We have the desired distribution. 


4. Use the inverse transformation method to generate a random variable having distribution 
function 

„, , x 2 + x 

F{x) = —-—, 0 < x < 1 

Solution: 


X 2 + X 


(x+lr-i 

X + - 
2 


X 


U 


2 U 


]/ 2U+ I 

\I 2U+ I 


1 

2 


(X, U G [0,1]) 
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By generating a random number, U, we have the desired distribution. 


U = rand ; 

X = sqrt ( 2 U + ^ 


5. Let X have a standard Cauchy distribution. 

Fx(x) = — arctan(x') H— 

7r 2 

Assuming you have U ~ Uniform^ 0,1), explain how to generate X. Then, use this result 
to produce 1000 samples of X and compute the sample mean. Repeat the experiment 100 
times. What do you observe and why? 

Solution: Using Inverse Transformation Method: 


U -= — arctan(A) 

2 7T 

7r ( U — = arctan(A) 


X = tan ( 7 t(U — -) 


Next, here is the MATLAB code: 


U = zeros( 1000,1); 
n = 100; 

average = zeros{n , 1); 
for i = 1 : n 
U = rand(1000,1); 

X = tan(pi * (U — 0.5)); 
average{i) = mean(X); 
end 

plotfaverage ) 


Cauchy distribution has no mean (Figure 12.6), or higher moments defined. 

6. (The Rejection Method) When we use the Inverse Transformation Method, we need a 
simple form of the cdf F(x) that allows direct computation of X = F 1_1 (t/). When 
F(x) doesn’t have a simple form but the pdf f(x) is available, random variables with 
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Figure 12.6: Cauchy Simulation 


density f(x) can be generated by the rejection method. Suppose you have a method 
for generating a random variable having density function g(x). Now, assume you want to 
generate a random variable having density function fix). Let c be a constant such that 

f(y) 


g(y) 


< c (for all y) 


Show that the following method generates a random variable with density function f{x). 

- Generate Y having density g. 

- Generate a random number U from Uniform (0,1). 

f(Y) 

-If U < ; set X = Y. Otherwise, return to step 1. 

Solution: The number of times N that the first two steps of the algorithm need to be called 
is itself a random variable and has a geometric distribution with “success” probability 

-OS) 

Thus, E{N ) = j-. Also, we can compute p: 

P\U < ^l\Y = y] = ^ 


cg(Y)' 


cg(y) 


p = 


f(y) 


I -oc eg {y) 

1 f°° 

- / f(y)dy 

C J —oc 


g(y)dy 


Therefore, E(N) = c 









16 CHAPTER 12. INTRODUCTION TO SIMULATION USING MATLABA. RAKHSHAN AND H. PISHRO-NIK 


Let F be the desired CDF (CDF of X). Now, we must show that the conditional distri¬ 
bution of Y given that U < ^g(Y) * s indeed F, i.e. P(Y < y\U < ) = F(y). Assume 

M = {U < cg(Y) I? K = {Y — v}- We know P{M) = p = \- Also, we can compute 

f(Y) P(U<l%X,Y<y) 

p(u<^\Y< y )= 

cg{Y) G(y) 


y P(U<^\Y = v<y) 


' —oo 
1 


G(y) 


g{v)dv 


_ - y f(v) 

G(y)J 

— OO cg(v) 

1 rv 


g{v)dv 


cG{y) J_ c 
F(y) 


f(v)dv 


Thus, 


cG(y) 

P(K\M) = P(M\K)P(K)/P{M) 
F{y) v G(y) 


cG(y) 

= F(y) 

7. Use the rejection method to generate a random variable having density function Beta(2,4). 
Hint: Assume g(x) = 1 for 0 < x < 1. 

Solution: 

f(x ) = 20x(l — x) 3 0 < x < 1 

g(x) = 1 0 < x < 1 

/(*) 




= 20x(l — x) 


We need to find the smallest constant c such that f(x)/g(x) < c. Differentiation of this 
quantity yields 


dx 

Thus, x = 


= 0 

1 
4 


Therefore, 


fix) 

135 
< - 

g{x) 

“ 64 

fix) 

256 

cgix) 

“ ~27 


Hence 
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n = 1 ; 

while(n == 1 ) 

U 1 = r and; 

U 2 = rand', 

if U2 <= 256/27 * C/l * (1 — C/1) A 3 

X = C/l; 

n = 0; 

end 

end 


8 . Use the rejection method to generate a random variable having the Gamma (§, 1) density 
function. Hint: Assume g(x) is the pdf of the Gamma (a = A = l). 

Solution: 


4 3 

f( x ) = 7r~i= x2 e X ,x>0 
oy 7r 

2 2a; 

f/(x) = -e“"5" x > 0 


fix) _ 

g(x) 3^ 

/O) 
fKU 


10 3 _ 3a: 

X 2 e 5 


dx 


= 0 


Hence, x = - 


c = 


/(s) 

c d( x ) 


10 / 5 

3 —3a: 

xae s 


51 I 


(I) 


e 2 


-3 

e 2 


We know how to generate an Exponential random variable. 


- Generate a random number U\ and set Y = — |logC/i. 

- Generate a random number C/ 2 . 

3 —3 y 

- If C /2 < ya 3 5 _ 3 , set A = y. Otherwise, execute the step 1 . 

(I) 2 ^ 


9. Use the rejection method to generate a standard normal random variable. Hint: Assume 
g(x) is the pdf of the exponential distribution with A = 1. 
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Solution: 


fix) = ~^=e 2 0 < x < oo 

g{x) = e~ x 0 < x < oo (Exponential density function with mean 1) 

f[x) l~2 A 
us, ——= \ - e 2 

g(x) V vr 

rpi , ■ ■ f{x) 

inus, x = 1 maximizes ——— 

9\ x ) 

Thus, c = \/— 


cg(x) 

- Generate Y, an exponential random variable with mean 1. 

- Generate a random number U. 

-(Y- 1) 2 

-If U < e 2 set X = Y. Otherwise, return to step 1. 

10. Use the rejection method to generate a Gammai 2,1) random variable conditional on its 
value being greater than 5, that is 


/(*) = 


xe~ x dx 


(x > 5) 


Hint: Assume g(x) be the density function of exponential distribution. 

Solution: Since Gamma(2,l) random variable has expected value 2, we use an exponential 
distribution with mean 2 that is conditioned to be greater than 5. 


/(*) = 


/ 5 °° xe(~ x ldx 


g(x) = 


J 5 

xeUO 

6e( _5 l 

1U-D 


x > 5 


x > 5 


ZM = * -(*=*) 

g{x) 3 

We obtain the maximum in x = 5 since is decreasing. Therefore, 

„ /(5) _ 5 


g( 5) 3 


- Generate a random number V. 
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- Y = 5-21og(y). 

- Generate a random number U. 

- If U < 2 \ set X = Y: otherwise return to step 1. 
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